{ "cells": [ { "cell_type": "markdown", "id": "c2f4fd6c", "metadata": { "pycharm": { "name": "#%% md\n" } }, "source": [ "# Implementing Shor's algorithm in Perceval" ] }, { "cell_type": "markdown", "id": "5162c9b8", "metadata": { "pycharm": { "name": "#%% md\n" } }, "source": [ "This notebook presents a simulation in Perceval of a 4-qubit 12-modes optical circuit performing Shor's algorithm, based on Alberto Politi, Jonathan C.F. Matthews, and Jeremy L. O'brien. \"Shor’s quantum factoring algorithm on a photonic chip.\" [Science](https://www.science.org/doi/10.1126/science.1173731) 325.5945 (2009): 1221-1221." ] }, { "cell_type": "markdown", "id": "f29e5abf", "metadata": { "pycharm": { "name": "#%% md\n" } }, "source": [ "## Shor's algorithm" ] }, { "cell_type": "markdown", "id": "5dd3a2a5", "metadata": { "pycharm": { "name": "#%% md\n" } }, "source": [ "The purpose of Shor's algorithm is to find nontrivial factors of a given number $N$ in polynomial time, yielding an near-exponential speedup compared to state of the art classical algortihms.\n", "\n", "The main routine of Shor's algorithm consists in finding the order $r$ of a number $a \\in \\mathbb{Z}_N$, i.e. the smallest integer $r$ such that $a^r = 1 \\pmod N$.\n", "\n", "If the order of a randoly chosen $a$ which is coprime with $N$ is even, then $(a^{r/2} - 1)(a^{r/2} + 1) = k N$. If none of these factors are multiples of $N$, then $\\gcd(N, a^{r/2} - 1)$ and $\\gcd(N, a^{r/2} + 1)$ are nontrivial factors of $N$." ] }, { "cell_type": "markdown", "id": "ff756302", "metadata": { "pycharm": { "name": "#%% md\n" } }, "source": [ "## Preliminaries" ] }, { "cell_type": "code", "execution_count": 13, "id": "8b013239", "metadata": { "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "import sys" ] }, { "cell_type": "code", "execution_count": 14, "id": "55fe7a18", "metadata": { "pycharm": { "name": "#%%\n" }, "scrolled": true }, "outputs": [], "source": [ "from IPython import display\n", "from collections import Counter\n", "from tabulate import tabulate\n", "from tqdm.auto import tqdm\n", "\n", "import sympy as sp\n", "import numpy as np\n", "from scipy.optimize import minimize\n", "\n", "import perceval as pcvl\n", "import quandelibc as qc\n", "import perceval.lib.phys as phys\n", "import perceval.lib.symb as symb" ] }, { "cell_type": "markdown", "id": "826377de", "metadata": { "pycharm": { "name": "#%% md\n" } }, "source": [ "### Path encoding functions\n", "\n", "The following functions allow for conversion between the qubit and Fock state representations." ] }, { "cell_type": "code", "execution_count": 3, "id": "cd2bbda8", "metadata": { "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "def toFockState(qubitState):\n", " # path encoding\n", " pe = {0:[1,0], 1:[0,1]}\n", " return [0] + pe[qubitState[0]] + pe[qubitState[2]] + [0, 0] + pe[qubitState[1]] + pe[qubitState[3]] + [0]" ] }, { "cell_type": "code", "execution_count": 4, "id": "55ee4ce8", "metadata": { "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "def toQubitState(fockState):\n", " # qubit modes\n", " x1 = [1, 2]\n", " f1 = [3, 4]\n", " x2 = [7, 8]\n", " f2 = [9, 10]\n", " # auxiliary modes\n", " am1 = [0, 5]\n", " am2 = [6, 11]\n", " \n", " # auxiliary modes\n", " for i in am1 + am2:\n", " if fockState[i]!=0:\n", " return None\n", " L=[]\n", " # qubit modes\n", " for q in [x1, x2, f1, f2]:\n", " if fockState[q[0]]+fockState[q[1]] != 1:\n", " return None\n", " else:\n", " L.append(fockState[q[1]])\n", " return L" ] }, { "cell_type": "code", "execution_count": 5, "id": "39b1c269", "metadata": { "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "def strState(state):\n", " return str(pcvl.BasicState(state))" ] }, { "cell_type": "markdown", "id": "016e2534", "metadata": { "pycharm": { "name": "#%% md\n" } }, "source": [ "## The circuit" ] }, { "cell_type": "markdown", "id": "01c60397", "metadata": { "pycharm": { "name": "#%% md\n" } }, "source": [ "### Quantum circuit\n", "\n", "The quantum circuit has been optimized after choosing parameters $N = 15$ and $a = 2$, and aims to calculate $r=4$.\n", "It features 5 qubits labelled $x_0, x_1, x_2$ and $f_1, f_2$. Qubits $x_i$ and $f_1$ are initially in state $|0\\rangle$, and $f_2$ in state $|1\\rangle$.\n", "In the non-optimised Shor algorithm, qubits $x_i$ encode a binary number representing a pre-image of the Modular Exponentiation Function (MEF) $x \\mapsto a^x \\pmod N$, while qubits $f_i$ hold the image obtained after applying the MEF to qubits $x_i$. Applying the MEF when qubits $x_i$ hold a superposition of different pre-images (obtained with H gates on qubits $x_i$) allows to efficiently compute the order $r$ of parameter $a$ modulo $N$.\n", "\n", "The circuit consists of $\\mathsf{H}$ gates being first applied to each $x_i$ qubit, followed by $\\mathsf{CNOT}$ gates applied on $x_1, f_1$ and $x_2, f_2$ pairs, where $x_i$ are control qubits; finally the inverse $\\mathsf{QFT}$ algorithm is applied on qubits $x_i$.\n", "\n", "$\\mathsf{CNOT}$ gates on $x_i, f_i$ pairs ($x_i$ being the control) are implemented using $\\mathsf{H}$ and $\\mathsf{CZ}$ gates: the $\\mathsf{CZ}$ gate is sandwiched between two applications of $\\mathsf{H}$ on $f_i$.\n", "\n", "The input state of the circuit after optimisation is $|0\\rangle_{x_0}|0\\rangle_{x_1}|0\\rangle_{x_2}|0\\rangle_{f_1}|1\\rangle_{f_2}$.\n", "\n", "The expected output state is then $\\frac{1}{2} |0\\rangle_{x_0} \\otimes \\left ( |0\\rangle_{x_1}|0\\rangle_{f_1} + |1\\rangle_{x_1}|1\\rangle_{f_1} \\right ) \\otimes \\left ( |0\\rangle_{x_2}|1\\rangle_{f_2} + |1\\rangle_{x_2}|0\\rangle_{f_2} \\right )$." ] }, { "cell_type": "markdown", "id": "b7b25cf5", "metadata": { "pycharm": { "name": "#%% md\n" } }, "source": [ "### Photonic circuit\n", "\n", "The optical circuit from the result by Politi et al features twelve modes (ordered from 0 to 11 from top to bottom).\n", "\n", "During the execution, qubit $x_0$ remains unentangled from the other qubits. It can therefore be removed from the optical implementation.\n", "\n", "The qubits $x_1, x_2, f_1, f_2$ are path encoded as modes $(1, 2)$, $(3, 4)$, $(7, 8)$, $(9, 10)$ respectively. The four remaining modes are used as auxiliary modes to implement the $\\mathsf{CZ}$ gates.\n", "\n", "With path encoding each $\\mathsf{H}$ gate in the quantum circuit is implemented with a beam splitter with reflectivity $R=1/2$ between the two pathes corresponding to the qubit. In our implementation in Perceval, phase shifters are added to properly tune the phase between each path.\n", "\n", "$\\mathsf{CZ}$ gates are implemented with three beam splitters with reflectivity $R=2/3$ acting on six modes: one inner BS creates interference between the two qubits, and two outer BS balance detection probability using auxiliary modes.\n", "This optical implementation succesfully yields the output state produced by a $\\mathsf{CZ}$ gate with probability 1/9; otherwise it creates a dummy state, which can be removed by post-selection.\n", "\n", "In the case $r=4$ the QFT can be performed classically and doesn't need to be implemented in the photonic circuit." ] }, { "cell_type": "markdown", "id": "3acb2273", "metadata": { "pycharm": { "name": "#%% md\n" } }, "source": [ "## In perceval" ] }, { "cell_type": "markdown", "id": "7aa60c25", "metadata": { "pycharm": { "name": "#%% md\n" } }, "source": [ "### Implementing the circuit" ] }, { "cell_type": "code", "execution_count": 6, "id": "2fe3a053", "metadata": { "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "circ = pcvl.Circuit(12)\n", "\n", "# qubit modes\n", "# for qubit states 0, 1\n", "x1 = [1, 2]\n", "f1 = [3, 4]\n", "x2 = [7, 8]\n", "f2 = [9, 10]\n", "# auxiliary modes\n", "am1 = [0, 5]\n", "am2 = [6, 11]\n", "\n", "\n", "# H gates\n", "for q in [x1, f1, x2, f2]:\n", " circ.add(q, symb.BS(R=1/2, phi=sp.pi/2))\n", " circ.add(q[1], symb.PS(phi=sp.pi))\n", "\n", "# CZ gates\n", "for x, f, am in [(x1, f1, am1), (x2, f2, am2)]:\n", " circ.add((am[0], x[0]), symb.BS(R=1/3)) # T = 1/3\n", " circ.add((x[1], f[0]), symb.BS(R=1/3))\n", " circ.add((f[1], am[1]), symb.BS(R=1/3))\n", "\n", "# H gates\n", "for q in [f1, f2]:\n", " circ.add(q, symb.BS(R=1/2, phi=sp.pi/2))\n", " circ.add(q[1], symb.PS(phi=sp.pi))" ] }, { "cell_type": "code", "execution_count": 7, "id": "034d5e45", "metadata": { "pycharm": { "name": "#%%\n" }, "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Optical circuit for Shor's algorithm\n" ] }, { "data": { "text/plain": "", "text/html": "\n0\n\n1\n\n2\n\n3\n\n4\n\n5\n\n6\n\n7\n\n8\n\n9\n\n10\n\n11\n\nΦ=pi/2\n\n\nΦ=pi\n\nΦ=pi/2\n\n\nΦ=pi\n\nΦ=pi/2\n\n\nΦ=pi\n\nΦ=pi/2\n\n\nΦ=pi\n\n\nR=1/3\n\n\nR=1/3\n\n\nR=1/3\n\n\nR=1/3\n\n\nR=1/3\n\n\nR=1/3\n\nΦ=pi/2\n\n\nΦ=pi\n\nΦ=pi/2\n\n\nΦ=pi\n\n\n\n\n\n\n\n\n\n\n\n0\n\n1\n\n2\n\n3\n\n4\n\n5\n\n6\n\n7\n\n8\n\n9\n\n10\n\n11" }, "metadata": {}, "output_type": "display_data" } ], "source": [ "print(\"Optical circuit for Shor's algorithm\")\n", "pcvl.pdisplay(circ)" ] }, { "cell_type": "code", "execution_count": 8, "id": "00d193de", "metadata": { "pycharm": { "name": "#%%\n" }, "scrolled": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The associated matrix\n" ] }, { "data": { "text/plain": "", "text/html": "$\\left[\\begin{array}{cccccccccccc}0.577350269189626 & 0.577350269189626 i & 0.577350269189626 i & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\0.816496580927726 i & 0.408248290463863 & 0.408248290463863 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\0 & 0.408248290463863 & -0.408248290463863 & 0.577350269189626 i & 0.577350269189626 i & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\0 & 0.408248290463863 i & - 0.408248290463863 i & 0.577350269189626 & 0 & 0.577350269189626 i & 0 & 0 & 0 & 0 & 0 & 0\\\\0 & 0.408248290463863 i & - 0.408248290463863 i & 0 & 0.577350269189626 & - 0.577350269189626 i & 0 & 0 & 0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0.577350269189626 i & - 0.577350269189626 i & 0.577350269189626 & 0 & 0 & 0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0 & 0 & 0 & 0.577350269189626 & 0.577350269189626 i & 0.577350269189626 i & 0 & 0 & 0\\\\0 & 0 & 0 & 0 & 0 & 0 & 0.816496580927726 i & 0.408248290463863 & 0.408248290463863 & 0 & 0 & 0\\\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0.408248290463863 & -0.408248290463863 & 0.577350269189626 i & 0.577350269189626 i & 0\\\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0.408248290463863 i & - 0.408248290463863 i & 0.577350269189626 & 0 & 0.577350269189626 i\\\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0.408248290463863 i & - 0.408248290463863 i & 0 & 0.577350269189626 & - 0.577350269189626 i\\\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0.577350269189626 i & - 0.577350269189626 i & 0.577350269189626\\end{array}\\right]$" }, "metadata": {}, "output_type": "display_data" } ], "source": [ "print(\"The associated matrix\")\n", "pcvl.pdisplay(circ.U)" ] }, { "cell_type": "markdown", "id": "7bdc227b", "metadata": { "pycharm": { "name": "#%% md\n" } }, "source": [ "### Input state" ] }, { "cell_type": "markdown", "id": "4e1ced97", "metadata": { "pycharm": { "name": "#%% md\n" } }, "source": [ "In the preliminaries we provided functions to map qubit states to the corresponding Fock states and vice-versa.\n", "\n", "A *computational basis qubit state* on qubits $x_1, f_1, x_2, f_2$ is represented with a list of 4 boolean values.\n", "\n", "A *Fock state* is represented with a list of twelve integer values.\n", "As described above, for Fock states, the modes are enumerated as follows:\n", "* mode pairs $(1,2), (3,4), (7,8), (9,10)$ respectively correspond to states $0,1$ for qubits $x_1, f_1, x_2, f_2$\n", "* modes $0,5,6,11$ are auxiliary modes used for CZ gates" ] }, { "cell_type": "markdown", "id": "51d9eca7", "metadata": { "pycharm": { "name": "#%% md\n" } }, "source": [ "The input state of the circuit is $|0\\rangle_{x_1}|0\\rangle_{x_2}|0\\rangle_{f_1}|1\\rangle_{f_2}$.\n", "With path encoding this corresponds to sending 4 photons in total in the optical circuit, in the input modes corresponding to the inital state of each qubit." ] }, { "cell_type": "code", "execution_count": 9, "id": "79f7ca45", "metadata": { "pycharm": { "name": "#%%\n" }, "scrolled": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Input qubit state: |0,0,0,1>\n", "Corresponding input Fock state: |0,1,0,1,0,0,0,1,0,0,1,0>\n" ] } ], "source": [ "qubit_istate = [0,0,0,1]\n", "istate = toFockState(qubit_istate)\n", "\n", "print(\"Input qubit state:\", strState(qubit_istate))\n", "print(\"Corresponding input Fock state:\", strState(istate))" ] }, { "cell_type": "markdown", "id": "d3f88c02", "metadata": { "pycharm": { "name": "#%% md\n" } }, "source": [ "## Simulation" ] }, { "cell_type": "code", "execution_count": 10, "id": "ca10ff5e", "metadata": { "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "backend = pcvl.BackendFactory().get_backend(\"Naive\")\n", "simulator = backend(circ)" ] }, { "cell_type": "markdown", "id": "e81a4f44", "metadata": { "pycharm": { "name": "#%% md\n" } }, "source": [ "### Computing the output state" ] }, { "cell_type": "markdown", "id": "c4e62d76", "metadata": { "pycharm": { "name": "#%% md\n" } }, "source": [ "Using perceval we compute the output state obtained with the optical circuit." ] }, { "cell_type": "markdown", "id": "2f3e6f17", "metadata": { "pycharm": { "name": "#%% md\n" } }, "source": [ "#### Amplitudes\n", "\n", "We first decompose the output state in the computational basis and plot the corresponding amplitudes." ] }, { "cell_type": "code", "execution_count": 11, "id": "0cd16853", "metadata": { "pycharm": { "name": "#%%\n" }, "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Output state amplitudes: (post-selected on qubit states, not renormalized)\n", "|x1,x2,f1,f2>\n", "|0,0,0,0> 0j\n", "|0,0,0,1> 0j\n", "|0,0,1,0> 0j\n", "|0,0,1,1> 0j\n", "|0,1,0,0> 0j\n", "|0,1,0,1> 0j\n", "|0,1,1,0> 0j\n", "|0,1,1,1> 0j\n", "|1,0,0,0> (-1.3877787807814457e-17+2.7755575615628914e-17j)\n", "|1,0,0,1> (-1.3877787807814457e-17+2.7755575615628914e-17j)\n", "|1,0,1,0> (-1.3877787807814457e-17+2.7755575615628914e-17j)\n", "|1,0,1,1> (-1.3877787807814457e-17+2.7755575615628914e-17j)\n", "|1,1,0,0> (-2.7755575615628914e-17+2.7755575615628914e-17j)\n", "|1,1,0,1> (-2.7755575615628914e-17+2.7755575615628914e-17j)\n", "|1,1,1,0> (-2.7755575615628914e-17+2.7755575615628914e-17j)\n", "|1,1,1,1> (-2.7755575615628914e-17+2.7755575615628914e-17j)\n" ] } ], "source": [ "output_qubit_states = [\n", " [x1,x2,f1,f2]\n", " for x1 in [0,1] for x2 in [0,1] for f1 in [0,1] for f2 in [0,1]\n", "]\n", "\n", "print(\"Output state amplitudes: (post-selected on qubit states, not renormalized)\")\n", "print(\"|x1,x2,f1,f2>\")\n", "for oqstate in output_qubit_states:\n", " ostate = toFockState(oqstate)\n", " a = simulator.probampli(pcvl.BasicState(pcvl.BasicState(istate)), pcvl.BasicState(ostate))\n", " print(strState(oqstate), a)" ] }, { "cell_type": "markdown", "id": "43e94f97", "metadata": { "pycharm": { "name": "#%% md\n" } }, "source": [ "The amplitudes obtained with perceval correspond to the expected output state\n", "$$\\frac{1}{2} \\left ( |0\\rangle_{x_1}|0\\rangle_{f_1} + |1\\rangle_{x_1}|1\\rangle_{f_1} \\right ) \\otimes \\left ( |0\\rangle_{x_2}|1\\rangle_{f_2} + |1\\rangle_{x_2}|0\\rangle_{f_2} \\right )$$\n", "up to numerical precision, without renormalization." ] }, { "cell_type": "markdown", "id": "300664ac", "metadata": { "pycharm": { "name": "#%% md\n" } }, "source": [ "#### Distribution\n", "\n", "We now compute the probabilities to obtain each outcome corresponding to measuring the expected output state in the computational basis." ] }, { "cell_type": "code", "execution_count": 12, "id": "2277ec12", "metadata": { "pycharm": { "name": "#%%\n" }, "scrolled": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Output state distribution: (post-selected on expected qubit states, not renormalized)\n", "|x1,x2,f1,f2>\n" ] }, { "data": { "text/plain": "", "text/html": "\n\n\n\n\n\n\n
|0,0,0,1> |0,1,0,0> |1,0,1,1> |1,1,1,0>
|0,0,0,1> 0 0 0 0
" }, "metadata": {}, "output_type": "display_data" } ], "source": [ "input_states = {\n", " pcvl.BasicState(pcvl.BasicState(istate)): strState(qubit_istate)\n", "}\n", "\n", "expected_output_states = {\n", " pcvl.BasicState(toFockState([x1,x2,x1,1-x2])): strState([x1,x2,x1,1-x2])\n", " for x1 in [0,1] for x2 in [0,1]\n", "}\n", "\n", "ca = pcvl.CircuitAnalyser(simulator, input_states, expected_output_states)\n", "ca.compute()\n", "\n", "print(\"Output state distribution: (post-selected on expected qubit states, not renormalized)\")\n", "print(\"|x1,x2,f1,f2>\")\n", "pcvl.pdisplay(ca)" ] }, { "cell_type": "markdown", "id": "5e464e2c", "metadata": { "pycharm": { "name": "#%% md\n" } }, "source": [ "The distribution computed with Perceval is uniform over each outcome, which corresponds to the expected distribution obtained in the paper when $x_0 = 0$." ] }, { "cell_type": "markdown", "source": [ "### Interpretation of the outcomes\n", "\n", "For each outcome we consider the values of qubits $x_2, x_1, x_0$ (where $x_0$ is 0) which represent a binary number between 0 and 7, here corresponding to 0, 4, 2 and 6 in the order of the table above.\n", "After sampling the circuit, obtainig outcomes 2 or 6 allows to successfully compute the order $r=4$.\n", "Obtaining outcome 0 is an expected failure of the quantum circuit, inherent to Shor's algorithm.\n", "Outcome 4 is an expected failure as well, as it only allows to compute the trivial factors 1 and 15.\n", "\n", "Since the distribution is uniform the circuit successfully yields a successful outcome with probability 1/2. This probability can be amplified exponentially close to $1$ by sampling the circuit multiple times." ], "metadata": { "collapsed": false, "pycharm": { "name": "#%% md\n" } } }, { "cell_type": "markdown", "source": [ "## Reference" ], "metadata": { "collapsed": false, "pycharm": { "name": "#%% md\n" } } }, { "cell_type": "markdown", "source": [ "Alberto Politi, Jonathan C.F. Matthews, and Jeremy L. O'brien. \"Shor’s quantum factoring algorithm on a photonic chip.\" [Science](https://www.science.org/doi/10.1126/science.1173731) 325.5945 (2009): 1221-1221." ], "metadata": { "collapsed": false, "pycharm": { "name": "#%% md\n" } } } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.4" } }, "nbformat": 4, "nbformat_minor": 5 }